• Descriptive point pattern analysis
• Kernel density maps
• Polygon extraction

We are going to use the tree species dataset from two weeks ago, more specifically the occurences of abies alba (silver fir):

library(sp)
library(raster)

abies <- read.csv("data/abiesalba.csv")

head(abies)

##          long      lat
## 1 -0.01377021 48.78623
## 2 -0.01423346 42.83379
## 3 -0.02027734 42.79684
## 4 -0.02039473 43.01494
## 5 -0.02171301 48.48641
## 6 -0.02336757 42.92377

abies_sp <- SpatialPoints(abies, CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))

abies_sp <- spTransform(abies_sp, CRS("+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"))

The centroid can easily be calculated by averaging the coordinates:

coords <- abies_sp@coords
coords <- as.data.frame(coords)

xmean <- mean(coords[, 1])
ymean <- mean(coords[, 2])

p <- ggplot(coords, aes(x = long, y = lat)) + 
  geom_point(alpha = 0.3) +
  geom_point(aes(x = xmean, y = ymean), col = "red", shape = 2)

We can also calculate the median of the coordinates:

xmedian <- median(coords[, 1])
ymedian <- median(coords[, 2])

p <- ggplot(coords, aes(x = long, y = lat)) + 
  geom_point(alpha = 0.3) +
  geom_point(aes(x = xmean, y = ymean), col = "red", shape = 2) +
  geom_point(aes(x = xmedian, y = ymedian), col = "blue", shape = 2)

We can calculate the average/minimum/maximum distance between points:

dist <- raster::pointDistance(abies_sp, longlat = FALSE)
dist[1:3, 1:3]

##          [,1]       [,2]       [,3]
## [1,]      0.0 661483.862 665579.331
## [2,] 661483.9      0.000   4123.066
## [3,] 665579.3   4123.066      0.000

dist <- dist[lower.tri(dist)]

c(min(dist), mean(dist), max(dist)) / 10000

## [1]   0.04996316  55.87651901 268.64010711

We can generate a grid from the points extent::

ext <- extent(abies_sp)
ext

## class       : Extent 
## xmin        : 3185500 
## xmax        : 5621500 
## ymin        : 1690500 
## ymax        : 3669500

ras <- raster(ext)
ras

## class       : RasterLayer 
## dimensions  : 10, 10, 100  (nrow, ncol, ncell)
## resolution  : 243600, 197900  (x, y)
## extent      : 3185500, 5621500, 1690500, 3669500  (xmin, xmax, ymin, ymax)
## coord. ref. : NA

And set the resolution to 1000*1000 meters while keeping the projection:

dim(ras) <- c(1000, 1000, 1)
ras

## class       : RasterLayer 
## dimensions  : 1000, 1000, 1e+06  (nrow, ncol, ncell)
## resolution  : 2436, 1979  (x, y)
## extent      : 3185500, 5621500, 1690500, 3669500  (xmin, xmax, ymin, ymax)
## coord. ref. : NA

projection(ras) <- projection(abies_sp)

Following we can use this raster to rasterize the points and count the number of points in each cell:

countgrid <- rasterize(abies_sp, ras, fun = "count", background = 0)

mapview(countgrid)

We can easily calculate the distances to cells with Abies. However, I first aggregate the raster to a loer resolution (otherwise it will take too long):

abies_ras <- aggregate(countgrid, fact = 10, fun = sum)

mapview(abies_ras)

Following I set all cells containing no Abies to NA and calculate the distance:

abies_ras[abies_ras == 0] <- NA

distance_abies <- distance(abies_ras)

mapview(distance_abies)

Kernel density estimates do not simply count the occurences, but 'smooth' by fitting a function (e.g., multivariate-normal) to the data. We can calculate Kernel densities as raster:

library(spatialEco)

kernel_ras <- sp.kde(abies_sp, bw = 10000, newdata = ras)

mapview(kernel_ras)

ggplot2 can also estimate kernels and nicely plot them atop of Google Maps!

library(ggmap)

map <- get_map(location = c(lon = mean(abies$long), lat = mean(abies$lat)), zoom = 5)

## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=47.077251,9.542565&zoom=5&size=640x640&scale=2&maptype=terrain&language=en-EN&sensor=false

p <- ggmap(map) + 
  geom_density2d(data = abies, aes(x = long, y = lat), size = 0.3) + 
  stat_density2d(data = abies, 
                 aes(x = long, y = lat, fill = ..level.., alpha = ..level..), size = 0.01, 
                 bins = 50, geom = "polygon") + 
  scale_fill_gradient(low = "green", high = "red") + 
  scale_alpha(range = c(0.1, 0.5), guide = FALSE)